Manuscript_CAPTAIN_analyses

Load libraries

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)
library(readr)
library(tidyr)
library(ggpubr)
library(betareg)
library(broom)
library(margins)
library(patchwork)
library(dplyr)
library(sf)
library(egg)
library(spatstat)
library(sf)
library(tidyverse)      # For data manipulation and visualization
library(rnaturalearthdata) # Additional natural earth data
library(biscale)        # For bivariate mapping
library(gridExtra)      # For arranging multiple plots
library(grid)           # For text elements in plots
library(colorspace)     # For color manipulation
library(knitr)
library(kableExtra)

Load the data

# Read the CAPTAIN2 EDGE2 RDS file
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 FUSE RDS file
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))

# Read the CAPTAIN2 IUCN RDS file
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
# Use the same raster file for consistency
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Read the protected range fractions RDS file
protected_fractions <- readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))

# Read the continental shark conservation metrics CSV file
shark_metrics <- read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))

# Load fishing data
load(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))

IUCN continental 0.1 budget

# Function to analyze priority data
analyze_priority_data <- function(data, index_name) {
  total_cells <- nrow(data)
  nonzero_cells <- sum(data$Priority > 0, na.rm = TRUE)
  zero_cells <- sum(data$Priority == 0, na.rm = TRUE)
  na_cells <- sum(is.na(data$Priority))
  
  # Get non-zero priorities for quantiles
  nonzero_priority <- data$Priority[data$Priority > 0]
  
  # Create summary data frame
  summary_df <- data.frame(
    Index = index_name,
    Total_Cells = total_cells,
    Nonzero_Cells = nonzero_cells,
    Nonzero_Percent = round(nonzero_cells/total_cells*100, 2),
    Zero_Cells = zero_cells,
    Zero_Percent = round(zero_cells/total_cells*100, 2),
    NA_Cells = na_cells,
    NA_Percent = round(na_cells/total_cells*100, 2),
    Min = min(data$Priority, na.rm = TRUE),
    Q25 = quantile(data$Priority, 0.25, na.rm = TRUE),
    Median = median(data$Priority, na.rm = TRUE),
    Mean = mean(data$Priority, na.rm = TRUE),
    Q75 = quantile(data$Priority, 0.75, na.rm = TRUE),
    Max = max(data$Priority, na.rm = TRUE),
    NonZero_Median = median(nonzero_priority, na.rm = TRUE),
    NonZero_Mean = mean(nonzero_priority, na.rm = TRUE)
  )
  
  return(summary_df)
}

# Analyze EDGE2 only
iucn_summary <- analyze_priority_data(CAPTAIN2_IUCN_data, "IUCN")

# Create a nicely formatted table
kable(iucn_summary,
      format = "html",
      digits = 2,
      col.names = c("Index", "Total Cells", "Non-zero Cells", "Non-zero %",
                    "Zero Cells", "Zero %", "NA Cells", "NA %",
                    "Min", "Q25", "Median", "Mean", "Q75", "Max",
                    "Non-zero Median", "Non-zero Mean"),
      caption = "Summary Statistics of IUCN Conservation Priority Values",
      row.names = FALSE) %>%  # Add this parameter
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 12) %>%
  add_header_above(c(" " = 1, "Cell Counts" = 7, "All Cells Statistics" = 6, "Non-zero Only" = 2)) %>%
  column_spec(1, bold = TRUE, width = "4em") %>%
  column_spec(2:8, width = "6em") %>%
  column_spec(9:16, width = "5em")

Summary Statistics of IUCN Conservation Priority Values

	Cell Counts							All Cells Statistics						Non-zero Only
Index	Total Cells	Non-zero Cells	Non-zero %	Zero Cells	Zero %	NA Cells	NA %	Min	Q25	Median	Mean	Q75	Max	Non-zero Median	Non-zero Mean
IUCN	232560	6031	2.59	226529	97.41	0	0	0	0	0	0.02	0	1	1	0.85

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 IUCN data based on PUID
CAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_IUCN_sf <- st_as_sf(
  CAPTAIN2_IUCN_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_IUCN <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_IUCN <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_IUCN <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_IUCN_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN

# Display the plot
print(CAPTAIN2_IUCN_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),
  plot = CAPTAIN2_IUCN_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

FUSE continental 0.1 budget

# Analyze EDGE2 only
fuse_summary <- analyze_priority_data(CAPTAIN2_FUSE_data, "FUSE")

# Create a nicely formatted table
kable(fuse_summary,
      format = "html",
      digits = 2,
      col.names = c("Index", "Total Cells", "Non-zero Cells", "Non-zero %",
                    "Zero Cells", "Zero %", "NA Cells", "NA %",
                    "Min", "Q25", "Median", "Mean", "Q75", "Max",
                    "Non-zero Median", "Non-zero Mean"),
      caption = "Summary Statistics of FUSE Conservation Priority Values",
      row.names = FALSE) %>%  # Add this parameter
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 12) %>%
  add_header_above(c(" " = 1, "Cell Counts" = 7, "All Cells Statistics" = 6, "Non-zero Only" = 2)) %>%
  column_spec(1, bold = TRUE, width = "4em") %>%
  column_spec(2:8, width = "6em") %>%
  column_spec(9:16, width = "5em")

Summary Statistics of FUSE Conservation Priority Values

	Cell Counts							All Cells Statistics						Non-zero Only
Index	Total Cells	Non-zero Cells	Non-zero %	Zero Cells	Zero %	NA Cells	NA %	Min	Q25	Median	Mean	Q75	Max	Non-zero Median	Non-zero Mean
FUSE	232560	7271	3.13	225289	96.87	0	0	0	0	0	0.02	0	1	1	0.75

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 FUSE data based on PUID
CAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_FUSE_sf <- st_as_sf(
  CAPTAIN2_FUSE_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_FUSE <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_FUSE <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_FUSE <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank())      # Add this line

# Create the plot
CAPTAIN2_FUSE_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE

# Display the plot
print(CAPTAIN2_FUSE_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),
  plot = CAPTAIN2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

EDGE2 continental 0.1 budget

# Analyze EDGE2 only
edge2_summary <- analyze_priority_data(CAPTAIN2_EDGE2_data, "EDGE2")

# Create a nicely formatted table
kable(edge2_summary,
      format = "html",
      digits = 2,
      col.names = c("Index", "Total Cells", "Non-zero Cells", "Non-zero %",
                    "Zero Cells", "Zero %", "NA Cells", "NA %",
                    "Min", "Q25", "Median", "Mean", "Q75", "Max",
                    "Non-zero Median", "Non-zero Mean"),
      caption = "Summary Statistics of EDGE2 Conservation Priority Values",
      row.names = FALSE) %>%  # Add this parameter
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 12) %>%
  add_header_above(c(" " = 1, "Cell Counts" = 7, "All Cells Statistics" = 6, "Non-zero Only" = 2)) %>%
  column_spec(1, bold = TRUE, width = "4em") %>%
  column_spec(2:8, width = "6em") %>%
  column_spec(9:16, width = "5em")

Summary Statistics of EDGE2 Conservation Priority Values

	Cell Counts							All Cells Statistics						Non-zero Only
Index	Total Cells	Non-zero Cells	Non-zero %	Zero Cells	Zero %	NA Cells	NA %	Min	Q25	Median	Mean	Q75	Max	Non-zero Median	Non-zero Mean
EDGE2	232560	5618	2.42	226942	97.58	0	0	0	0	0	0.02	0	1	1	0.96

# Load one of your input raster files to extract the correct grid structure
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 data based on PUID
CAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_EDGE2_sf <- st_as_sf(
  CAPTAIN2_EDGE2_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_EDGE2 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_EDGE2 <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_EDGE2 <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Create the plot
CAPTAIN2_EDGE2_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the plot
print(CAPTAIN2_EDGE2_plot)

# Save the plot
ggsave(
  filename = here::here("Outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),
  plot = CAPTAIN2_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Difference maps

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Join all three datasets with coordinates
IUCN_with_coords <- CAPTAIN2_IUCN_data %>%
  dplyr::select(PUID, IUCN = Priority) %>%
  left_join(coords, by = "PUID")

EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>%
  dplyr::select(PUID, EDGE2 = Priority) %>%
  left_join(coords, by = "PUID") 

FUSE_with_coords <- CAPTAIN2_FUSE_data %>%
  dplyr::select(PUID, FUSE = Priority) %>%
  left_join(coords, by = "PUID")

# Combine all datasets
all_indices <- coords %>%
  left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by = "PUID")

# Calculate differences
all_indices <- all_indices %>%
  mutate(
    # Replace NA with 0 for calculation purposes
    IUCN = ifelse(is.na(IUCN), 0, IUCN),
    EDGE2 = ifelse(is.na(EDGE2), 0, EDGE2),
    FUSE = ifelse(is.na(FUSE), 0, FUSE),
    
    # Calculate differences
    IUCN_minus_FUSE = IUCN - FUSE,
    IUCN_minus_EDGE2 = IUCN - EDGE2,
    EDGE2_minus_FUSE = EDGE2 - FUSE
  )

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Filter to non-zero differences for each comparison to reduce plot size
# IUCN - FUSE
IUCN_FUSE_diff <- all_indices %>%
  filter(IUCN_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# IUCN - EDGE2
IUCN_EDGE2_diff <- all_indices %>%
  filter(IUCN_minus_EDGE2 != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# EDGE2 - FUSE
EDGE2_FUSE_diff <- all_indices %>%
  filter(EDGE2_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create a diverging color palette for difference maps
# Blue for negative (first index lower), white for zero, red for positive (first index higher)
diff_colors <- c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", 
                "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Display all plots
print(IUCN_FUSE_plot)

print(IUCN_EDGE2_plot)

print(EDGE2_FUSE_plot)

# Save all plots
ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference_2.png"),
  plot = IUCN_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference_2.png"),
  plot = IUCN_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference_2.png"),
  plot = EDGE2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Species level priorities

# Rename column in shark_metrics to match better
shark_metrics <- shark_metrics %>%
  rename(Species = `Species name`)

# Order shark_metrics alphabetically by Species name
shark_metrics <- shark_metrics %>%
  arrange(Species)

# Add an original order ID to protected_fractions to maintain its original order
protected_fractions$original_order <- 1:nrow(protected_fractions)

# Add row number as IDs to both datasets
protected_fractions$protected_ID <- 1:nrow(protected_fractions)
shark_metrics$species_ID <- 1:nrow(shark_metrics)

# Check if the datasets have the same number of rows
if(nrow(protected_fractions) == nrow(shark_metrics)) {
  # Create index mapping
  indices <- data.frame(
    protected_ID = 1:nrow(protected_fractions),
    species_ID = 1:nrow(shark_metrics)
  )
  
  # Join protected_fractions with indices
  protected_with_indices <- protected_fractions %>%
    left_join(indices, by = "protected_ID")
  
  # Join shark_metrics with indices
  shark_with_indices <- shark_metrics %>%
    left_join(indices, by = "species_ID")
  
  # Now join the datasets
  combined_data <- protected_with_indices %>%
    inner_join(
      shark_with_indices,
      by = c("species_ID", "protected_ID"),
      suffix = c("_captain", "_original")
    ) %>%
    arrange(original_order)
  
  cat("Successfully joined datasets with", nrow(combined_data), "species\n")
  
  # Define IUCN categories and order
  iucn_order <- c("LC", "NT", "VU", "EN", "CR")
  
  iucn_colors <- c(
    "LC" = "#50C878",
    "NT" = "#FFFF00",
    "VU" = "#FFA500",
    "EN" = "#FF8C00",
    "CR" = "#FF0000"
  )
  
  # Calculate sample sizes for each group
  iucn_n <- combined_data %>%
    group_by(IUCN_original) %>%
    summarise(n = n()) %>%
    arrange(IUCN_original)
  
  fuse_n <- combined_data %>%
    group_by(FUSE_original) %>%
    summarise(n = n()) %>%
    arrange(FUSE_original)
  
  edge2_n <- combined_data %>%
    group_by(EDGE2_original) %>%
    summarise(n = n()) %>%
    arrange(EDGE2_original)
  
  # Create x-axis labels with sample sizes
  iucn_labels <- paste0(iucn_order, "\n(n=", iucn_n$n, ")")
  fuse_labels <- paste0(1:5, "\n(n=", fuse_n$n, ")")
  edge2_labels <- paste0(1:5, "\n(n=", edge2_n$n, ")")
  
  # 1. IUCN plot
  iucn_plot <- combined_data %>%
    mutate(
      IUCN_status = factor(IUCN_original, levels = 1:5, labels = iucn_order),
      protection_percentage = IUCN_captain * 100
    ) %>%
    ggplot(aes(x = IUCN_status, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "IUCN Red List threat status",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(limits = iucn_order, labels = iucn_labels) +
    scale_fill_manual(values = iucn_colors) +
    scale_color_manual(values = iucn_colors) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 2. FUSE plot
  fuse_plot <- combined_data %>%
    mutate(
      FUSE_category = factor(FUSE_original),
      protection_percentage = FUSE_captain * 100
    ) %>%
    ggplot(aes(x = FUSE_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "FUSE Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = fuse_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # 3. EDGE2 plot
  edge2_plot <- combined_data %>%
    mutate(
      EDGE2_category = factor(EDGE2_original),
      protection_percentage = EDGE2_captain * 100
    ) %>%
    ggplot(aes(x = EDGE2_category, y = protection_percentage)) +
    geom_jitter(width = 0.1,
                size = 0.4,
                alpha = 0.5,
                color = "darkgray") +
    stat_summary(fun = mean,
                 geom = "point",
                 size = 3,
                 color = "black") +
    stat_summary(fun.data = function(x) {
      mean_val <- mean(x, na.rm = TRUE)
      sd_val <- sd(x, na.rm = TRUE)
      return(data.frame(y = mean_val,
                        ymin = max(0, mean_val - sd_val),
                        ymax = min(100, mean_val + sd_val)))
    },
    geom = "errorbar",
    width = 0.1,
    color = "black") +
    labs(
      x = "EDGE2 Score",
      y = "Range protected (%)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
      legend.position = "none",
      panel.grid.major.x = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      axis.title.x = element_text(size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x = element_text(size = 14),
      axis.text.y = element_text(size = 14)
    ) +
    scale_x_discrete(labels = edge2_labels) +
    scale_y_continuous(limits = c(0, 100),
                       breaks = seq(0, 100, 25))
  
  # Combine the three plots into a grid
  combined_plots <- ggpubr::ggarrange(
    iucn_plot,
    fuse_plot,
    edge2_plot,
    ncol = 3,
    nrow = 1,
    labels = c("(A)", "(B)", "(C)"),
    font.label = list(size = 14,
                      face = "bold",
                      hjust = -2,
                      vjust = 1.5)
  )
  
  #print(combined_plots)
  
  # Save combined plot
  ggsave(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"),
         combined_plots, width = 20, height = 6, dpi = 150, bg = "white")

  #======================
  # Beta Regression Analysis
  #======================
  
  # Prepare data for beta regression
  model_data <- combined_data %>%
    mutate(
      IUCN_status = IUCN_original,
      FUSE_score = FUSE_original,
      EDGE2_score = EDGE2_original,
      IUCN_protection = IUCN_captain * 100,
      FUSE_protection = FUSE_captain * 100,
      EDGE2_protection = EDGE2_captain * 100,
      # Convert to (0,1) scale and handle boundaries
      n_obs = n(),
      IUCN_protection_beta = (IUCN_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      FUSE_protection_beta = (FUSE_protection/100 * (n_obs - 1) + 0.5) / n_obs,
      EDGE2_protection_beta = (EDGE2_protection/100 * (n_obs - 1) + 0.5) / n_obs
    )
  
  # Fit beta regression models
  iucn_beta <- betareg(IUCN_protection_beta ~ IUCN_status, data = model_data)
  fuse_beta <- betareg(FUSE_protection_beta ~ FUSE_score, data = model_data)
  edge2_beta <- betareg(EDGE2_protection_beta ~ EDGE2_score, data = model_data)
  
  # View summaries
  cat("\n=== IUCN Beta Regression ===\n")
  print(summary(iucn_beta))
  
  cat("\n=== FUSE Beta Regression ===\n")
  print(summary(fuse_beta))
  
  cat("\n=== EDGE2 Beta Regression ===\n")
  print(summary(edge2_beta))
  
  # Get tidy summaries
  iucn_tidy <- tidy(iucn_beta)
  fuse_tidy <- tidy(fuse_beta)
  edge2_tidy <- tidy(edge2_beta)
  
  # Get pseudo R²
  iucn_r2 <- cor(predict(iucn_beta), model_data$IUCN_protection_beta)^2
  fuse_r2 <- cor(predict(fuse_beta), model_data$FUSE_protection_beta)^2
  edge2_r2 <- cor(predict(edge2_beta), model_data$EDGE2_protection_beta)^2
  
  # Calculate marginal effects (average change in percentage per unit increase)
  iucn_margins <- margins(iucn_beta)
  fuse_margins <- margins(fuse_beta)
  edge2_margins <- margins(edge2_beta)
  
  # Extract AME and SE
  iucn_ame_summary <- summary(iucn_margins)
  fuse_ame_summary <- summary(fuse_margins)
  edge2_ame_summary <- summary(edge2_margins)
  
  # Convert to percentage (margins gives proportions)
  iucn_ame <- iucn_ame_summary$AME * 100
  iucn_ame_se <- iucn_ame_summary$SE * 100
  
  fuse_ame <- fuse_ame_summary$AME * 100
  fuse_ame_se <- fuse_ame_summary$SE * 100
  
  edge2_ame <- edge2_ame_summary$AME * 100
  edge2_ame_se <- edge2_ame_summary$SE * 100
  
  cat("\n=== Average Marginal Effects (as percentages) ===\n")
  cat(paste("IUCN AME:", round(iucn_ame, 2), "%, SE:", round(iucn_ame_se, 2), "%\n"))
  cat(paste("FUSE AME:", round(fuse_ame, 2), "%, SE:", round(fuse_ame_se, 2), "%\n"))
  cat(paste("EDGE2 AME:", round(edge2_ame, 2), "%, SE:", round(edge2_ame_se, 2), "%\n"))
  
  # Calculate predicted values and differences
  # IUCN predictions
  iucn_pred <- predict(iucn_beta,
                       newdata = data.frame(IUCN_status = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by IUCN category:\n")
  print(data.frame(Category = 1:5,
                   IUCN_Category = c("LC", "NT", "VU", "EN", "CR"),
                   Predicted_Protection = round(iucn_pred, 2)))
  
  iucn_avg_increase <- mean(diff(iucn_pred))
  cat(paste("Average increase per IUCN category:", round(iucn_avg_increase, 2), "%\n"))
  
  # FUSE predictions
  fuse_pred <- predict(fuse_beta,
                       newdata = data.frame(FUSE_score = 1:5),
                       type = "response") * 100
  
  cat("\nPredicted protection by FUSE score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(fuse_pred, 2)))
  
  fuse_avg_increase <- mean(diff(fuse_pred))
  cat(paste("Average increase per FUSE unit:", round(fuse_avg_increase, 2), "%\n"))
  
  # EDGE2 predictions
  edge2_pred <- predict(edge2_beta,
                        newdata = data.frame(EDGE2_score = 1:5),
                        type = "response") * 100
  
  cat("\nPredicted protection by EDGE2 score:\n")
  print(data.frame(Score = 1:5,
                   Predicted_Protection = round(edge2_pred, 2)))
  
  edge2_avg_increase <- mean(diff(edge2_pred))
  cat(paste("Average increase per EDGE2 unit:", round(edge2_avg_increase, 2), "%\n"))
  
  # Create comprehensive summary table with AME
  comprehensive_summary <- data.frame(
    Approach = c("IUCN", "FUSE", "EDGE2"),
    Coefficient_logit = c(
      iucn_tidy$estimate[2],
      fuse_tidy$estimate[2],
      edge2_tidy$estimate[2]
    ),
    SE = c(
      iucn_tidy$std.error[2],
      fuse_tidy$std.error[2],
      edge2_tidy$std.error[2]
    ),
    z_value = c(
      iucn_tidy$statistic[2],
      fuse_tidy$statistic[2],
      edge2_tidy$statistic[2]
    ),
    P_value = c(
      iucn_tidy$p.value[2],
      fuse_tidy$p.value[2],
      edge2_tidy$p.value[2]
    ),
    AME_percent = c(
      iucn_ame,
      fuse_ame,
      edge2_ame
    ),
    AME_SE = c(
      iucn_ame_se,
      fuse_ame_se,
      edge2_ame_se
    ),
    Pseudo_R2 = c(iucn_r2, fuse_r2, edge2_r2)
  )
  
  cat("\n=== Comprehensive Beta Regression Summary ===\n")
  print(comprehensive_summary)
  
  # Save results
  write.csv(comprehensive_summary,
            here::here("Outputs", "beta_regression_summary.csv"),
            row.names = FALSE)
  
} else {
  cat("ERROR: Datasets have different number of rows.\n")
  cat("Protected fractions:", nrow(protected_fractions), "rows\n")
  cat("Shark metrics:", nrow(shark_metrics), "rows\n")
}

Successfully joined datasets with 1000 species

=== IUCN Beta Regression ===

Call:
betareg(formula = IUCN_protection_beta ~ IUCN_status, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5313 -0.5004 -0.0020  0.5608  2.0202 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.05978    0.07511  -0.796    0.426    
IUCN_status  0.11768    0.02988   3.938 8.21e-05 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.91473    0.03196   28.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 379.5 on 3 Df
Pseudo R-squared: 0.01504
Number of iterations: 7 (BFGS) + 1 (Fisher scoring) 

=== FUSE Beta Regression ===

Call:
betareg(formula = FUSE_protection_beta ~ FUSE_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.3375 -0.4546  0.0051  0.5057  2.3049 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.19664    0.08966  -2.193   0.0283 *
FUSE_score   0.14261    0.06830   2.088   0.0368 *

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  1.04726    0.03673   28.52   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 211.9 on 3 Df
Pseudo R-squared: 0.004321
Number of iterations: 7 (BFGS) + 2 (Fisher scoring) 

=== EDGE2 Beta Regression ===

Call:
betareg(formula = EDGE2_protection_beta ~ EDGE2_score, data = model_data)

Quantile residuals:
    Min      1Q  Median      3Q     Max 
-2.5490 -0.4954 -0.0052  0.4548  2.1662 

Coefficients (mean model with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.68964    0.09615  -7.173 7.34e-13 ***
EDGE2_score  0.46406    0.07688   6.036 1.58e-09 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)  0.84174    0.02926   28.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 506.7 on 3 Df
Pseudo R-squared: 0.04761
Number of iterations: 13 (BFGS) + 1 (Fisher scoring) 

=== Average Marginal Effects (as percentages) ===
IUCN AME: 2.9 %, SE: 0.73 %
FUSE AME: 3.56 %, SE: 1.7 %
EDGE2 AME: 11.39 %, SE: 1.85 %

Predicted protection by IUCN category:
  Category IUCN_Category Predicted_Protection
1        1            LC                51.45
2        2            NT                54.38
3        3            VU                57.28
4        4            EN                60.13
5        5            CR                62.92
Average increase per IUCN category: 2.87 %

Predicted protection by FUSE score:
  Score Predicted_Protection
1     1                48.65
2     2                52.21
3     3                55.75
4     4                59.24
5     5                62.63
Average increase per FUSE unit: 3.5 %

Predicted protection by EDGE2 score:
  Score Predicted_Protection
1     1                44.38
2     2                55.93
3     3                66.88
4     4                76.25
5     5                83.63
Average increase per EDGE2 unit: 9.81 %

=== Comprehensive Beta Regression Summary ===
            Approach Coefficient_logit         SE  z_value      P_value
IUCN_status     IUCN         0.1176756 0.02987957 3.938332 8.205012e-05
FUSE_score      FUSE         0.1426108 0.06830262 2.087925 3.680458e-02
EDGE2_score    EDGE2         0.4640609 0.07687772 6.036350 1.576390e-09
            AME_percent    AME_SE  Pseudo_R2
IUCN_status    2.898230 0.7271499 0.01750165
FUSE_score     3.558461 1.6986270 0.00239858
EDGE2_score   11.386873 1.8514148 0.02982912

  # Display the saved image in the Quarto document
knitr::include_graphics(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"))

Manuscript maps

Individual maps

CAPTAIN2_EDGE2_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority EDGE2",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

CAPTAIN2_FUSE_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority FUSE",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE

CAPTAIN2_IUCN_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority IUCN",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid1_maps_continental_2ndrun_new.png"), grid1, width = 8, height = 15, dpi = 300)

Difference maps

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority (EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "A")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "B")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Save the grids if needed
ggsave(here::here("Outputs", "grid3_maps_continental_2ndrun.png"), grid2, width = 16, height = 15, dpi = 300)

Combined individual and difference maps

library(patchwork)

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Second grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "D")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "E")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "F")

# Create each grid (3 rows, 1 column)
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Combine the two grids side by side (2 columns)
combined_grid <- grid1 | grid2

# If you want to save it
ggsave(
  filename = here::here("Outputs", "Fig.3_combined_priority_maps.png"),
  plot = combined_grid,
  width = 10,    # Adjust width as needed for two columns
  height = 12,   # Adjust height as needed for three rows
  dpi = 300,
  bg = "white"
)

# Display the saved image in the Quarto document
knitr::include_graphics(here::here("Outputs", "Fig.3_combined_priority_maps.png"))

Spatial distribution analysis

# Filter for high priority cells (>0.9) and calculate percentages
# Calculate total coastal cells and high-priority cells for each index
# IUCN
CAPTAIN2_IUCN_high <- CAPTAIN2_IUCN_sf %>% 
  filter(Priority > 0.9)
iucn_total_cells <- nrow(CAPTAIN2_IUCN_sf)
iucn_high_cells <- nrow(CAPTAIN2_IUCN_high)
iucn_percent <- round((iucn_high_cells / iucn_total_cells) * 100, 1)

# FUSE
CAPTAIN2_FUSE_high <- CAPTAIN2_FUSE_sf %>% 
  filter(Priority > 0.9)
fuse_total_cells <- nrow(CAPTAIN2_FUSE_sf)
fuse_high_cells <- nrow(CAPTAIN2_FUSE_high)
fuse_percent <- round((fuse_high_cells / fuse_total_cells) * 100, 1)

# EDGE2
CAPTAIN2_EDGE2_high <- CAPTAIN2_EDGE2_sf %>% 
  filter(Priority > 0.9)
edge2_total_cells <- nrow(CAPTAIN2_EDGE2_sf)
edge2_high_cells <- nrow(CAPTAIN2_EDGE2_high)
edge2_percent <- round((edge2_high_cells / edge2_total_cells) * 100, 1)


#Mean nearest neighbor for high priority grid cells:

# Function to calculate Mean Nearest Neighbor Distance
calculate_mnn <- function(sf_high_priority) {
  # Get centroids and convert to coordinates
  coords <- st_coordinates(st_centroid(sf_high_priority))
  
  # Create a bounding box for the window
  bbox <- st_bbox(sf_high_priority)
  
  # Create observation window
  window <- owin(xrange = c(bbox["xmin"], bbox["xmax"]),
                 yrange = c(bbox["ymin"], bbox["ymax"]))
  
  # Create point pattern object
  ppp_obj <- ppp(coords[,1], coords[,2], window = window)
  
  # Calculate mean nearest neighbor distance
  mnn <- mean(nndist(ppp_obj))
  
  # Calculate standard deviation for context
  mnn_sd <- sd(nndist(ppp_obj))
  
  return(list(mean = mnn, sd = mnn_sd))
}

# Calculate for each index
iucn_mnn <- calculate_mnn(CAPTAIN2_IUCN_high)
fuse_mnn <- calculate_mnn(CAPTAIN2_FUSE_high)
edge2_mnn <- calculate_mnn(CAPTAIN2_EDGE2_high)

# Create summary data frame
mnn_summary <- data.frame(
  Index = c("IUCN", "FUSE", "EDGE2"),
  N_cells = c(iucn_high_cells, fuse_high_cells, edge2_high_cells),
  Percent = c(iucn_percent, fuse_percent, edge2_percent),
  MNN_mean = c(iucn_mnn$mean, fuse_mnn$mean, edge2_mnn$mean),
  MNN_sd = c(iucn_mnn$sd, fuse_mnn$sd, edge2_mnn$sd)
)

print(mnn_summary)

  Index N_cells Percent MNN_mean    MNN_sd
1  IUCN    4734    78.5 51143.37 17182.874
2  FUSE    4396    60.5 49594.94 26765.936
3 EDGE2    5328    94.8 48211.63  8975.762

Congruence maps

# Method 1: If you have separate dataframes for each index
# Assuming your data has columns: PUID, Priority, and geometry

# Fxirst, identify high priority areas (>0.9) for each index
high_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority > 0.9, ]
high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority > 0.9, ]
high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority > 0.9, ]

# Find congruent areas (present in all three)
# Method using PUID (assuming you have PUID column)
congruent_PUIDs <- intersect(intersect(high_priority_EDGE2$PUID, 
                                      high_priority_FUSE$PUID), 
                            high_priority_IUCN$PUID)

# Extract congruent areas
congruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]

# Method 2: If you need to merge dataframes first
# Create a combined dataframe with all three priority scores
# First, extract just the data (without geometry) from the other sf objects
FUSE_data <- st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])
IUCN_data <- st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])

# Merge with the EDGE2 sf object (keeping geometry)
combined_priorities <- merge(CAPTAIN2_EDGE2_sf, FUSE_data, 
                            by = "PUID", suffixes = c("_EDGE2", "_FUSE"))
combined_priorities <- merge(combined_priorities, IUCN_data, 
                            by = "PUID")
names(combined_priorities)[names(combined_priorities) == "Priority"] <- "Priority_IUCN"

# Identify congruent high priority areas
congruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 > 0.9 & 
                                         combined_priorities$Priority_FUSE > 0.9 & 
                                         combined_priorities$Priority_IUCN > 0.9, ]

# Create a map showing only congruent areas
congruent_map <- ggplot() +
  geom_sf(data = congruent_areas, aes(fill = "Congruent High Priority"), 
          color = "red", size = 0.8, alpha = 0.8) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_fill_manual(values = c("Congruent high priority" = "red"),
                    name = "Priority areas") +
  labs(title = "Congruent high priority areas",
       subtitle = "Areas with priority > 0.9 in IUCN, FUSE and EDGE2 dimensions",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the map
print(congruent_map)

ggsave(here::here("Outputs", "congruent_high_prioritiy_areas_09.png"), 
       congruent_map, 
       width = 10, height = 8, dpi = 300, bg = "white")

# Create detailed congruence categories
# Define logical conditions for each index
EDGE2_high <- combined_priorities$Priority_EDGE2 > 0.9
FUSE_high <- combined_priorities$Priority_FUSE > 0.9
IUCN_high <- combined_priorities$Priority_IUCN > 0.9

# Create detailed congruence categories
combined_priorities$congruence_category <- case_when(
  EDGE2_high & FUSE_high & IUCN_high ~ "All three indices",
  !EDGE2_high & FUSE_high & IUCN_high ~ "IUCN + FUSE",      # Changed from "FUSE + IUCN"
  EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN + EDGE2",     # Changed from "EDGE2 + IUCN"
  EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE + EDGE2",     # Changed from "EDGE2 + FUSE"
  !EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN only",
  !EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE only",
  EDGE2_high & !FUSE_high & !IUCN_high ~ "EDGE2 only",
  TRUE ~ "No high priority"
)

# Convert to factor with desired order
combined_priorities$congruence_category <- factor(
  combined_priorities$congruence_category,
  levels = c("All three indices",
             "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed order
             "IUCN only", "FUSE only", "EDGE2 only")
)

# Filter data to only include high priority areas
high_priority_data <- combined_priorities[combined_priorities$congruence_category %in%
                                          c("All three indices", 
                                            "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2",  # Changed names
                                            "IUCN only", "FUSE only", "EDGE2 only"), ]

# If simple plot fails, let's check the data more thoroughly
if (nrow(high_priority_data) == 0) {
  cat("No high priority data found! Checking original data...\n")
  cat("EDGE2 > 0.9:", sum(combined_priorities$Priority_EDGE2 > 0.9, na.rm = TRUE), "\n")
  cat("FUSE > 0.9:", sum(combined_priorities$Priority_FUSE > 0.9, na.rm = TRUE), "\n") 
  cat("IUCN > 0.9:", sum(combined_priorities$Priority_IUCN > 0.9, na.rm = TRUE), "\n")
}

# Extract coordinates for the congruence plot
if (nrow(high_priority_data) > 0) {
  coords <- st_coordinates(st_centroid(high_priority_data))
  plot_data <- data.frame(
    x = coords[,1],
    y = coords[,2], 
    category = high_priority_data$congruence_category
  )
  
  # Remove any rows with NA category
  plot_data <- plot_data[!is.na(plot_data$category), ]
  
  #cat("Creating enhanced congruence plot...\n")
  #cat("Plot data dimensions:", nrow(plot_data), "points\n")
  
  # Create the enhanced congruence plot
  congruence_map <- ggplot(plot_data, aes(x = x, y = y, color = category)) +
    geom_point(size = 1.2, alpha = 0.85, stroke = 0) +  # Larger points, no stroke to reduce overlap
    scale_color_manual(
      values = c("All three indices" = "#8B0000",      # Dark red
                "EDGE2 + FUSE" = "#FF8C00",            # Dark orange  
                "EDGE2 + IUCN" = "#DC143C",            # Crimson
                "FUSE + IUCN" = "#FFD700",             # Gold
                "EDGE2 only" = "#4169E1",              # Royal blue
                "FUSE only" = "#32CD32",               # Lime green
                "IUCN only" = "#9370DB"),              # Medium purple
      name = "High Priority\nCongruence",
      guide = guide_legend(
        override.aes = list(size = 4, alpha = 1),  # Larger, more opaque legend points
        title.position = "top",
        title.hjust = 0.5,
        ncol = 4  # Two columns for legend to save space
      )) +
    labs(
      title = "Global Conservation Priority Congruence Analysis",
      subtitle = "High priority areas (>0.9) showing agreement patterns across EDGE2, FUSE, and IUCN indices",
      x = "Longitude", 
      y = "Latitude",
      caption = paste("Total high priority areas:", nrow(plot_data))
    ) +
    theme_minimal() +
    theme(
      # Plot aesthetics
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "#f8f9fa", color = NA),
      panel.grid.major = element_line(color = "white", size = 0.5, linetype = "solid"),
      panel.grid.minor = element_line(color = "white", size = 0.25, linetype = "solid"),
      
      # Text styling
      plot.title = element_text(size = 16, hjust = 0.5, face = "bold", 
                               margin = margin(b = 10)),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray30",
                                  margin = margin(b = 20)),
      plot.caption = element_text(size = 10, color = "gray50", hjust = 1),
      
      # Axis styling
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(size = 10, color = "gray30"),
      axis.ticks = element_line(color = "gray50", size = 0.5),
      
      # Legend styling
      legend.position = "bottom",
      legend.background = element_rect(fill = "white", color = "gray80", size = 0.5),
      legend.margin = margin(t = 15, r = 10, b = 10, l = 10),
      legend.title = element_text(size = 12, face = "bold"),
      legend.text = element_text(size = 10),
      legend.key.size = unit(0.8, "cm"),
      
      # Panel border
      panel.border = element_rect(color = "gray70", fill = NA, size = 0.5)
    ) +
    # Add coordinate system and aspect ratio
    coord_fixed(ratio = 1.3) +  # Adjust ratio for better map appearance
    # Add subtle borders around plot area
    scale_x_continuous(expand = c(0.02, 0.02)) +
    scale_y_continuous(expand = c(0.02, 0.02))
  
  # Print the enhanced plot
  #print(congruence_map)
  #cat("Enhanced congruence plot created successfully!\n")
  
  # Print summary by category
  cat("\nBreakdown by congruence category:\n")
  category_counts <- table(plot_data$category)
  category_percentages <- round(prop.table(category_counts) * 100, 1)
  
  for(i in 1:length(category_counts)) {
    cat(sprintf("%-20s: %4d points (%4.2f%%)\n", 
                names(category_counts)[i], 
                category_counts[i], 
                category_percentages[i]))
  }
  
} else {
  cat("No valid high priority data to plot\n")
}

Breakdown by congruence category:
All three indices   : 1459 points (46.50%)
IUCN + FUSE         :   75 points (2.40%)
IUCN + EDGE2        : 1130 points (36.00%)
FUSE + EDGE2        :  275 points (8.80%)
IUCN only           :   29 points (0.90%)
FUSE only           :   17 points (0.50%)
EDGE2 only          :  154 points (4.90%)

# Print summary of congruence patterns
#cat("\nCongruence Pattern Summary:\n")
congruence_summary <- table(combined_priorities$congruence_category, useNA = "ifany")
#print(congruence_summary)

# Calculate percentages
total_high_priority <- sum(congruence_summary[names(congruence_summary) != "No high priority"], na.rm = TRUE)
congruence_percentages <- round(congruence_summary / total_high_priority * 100, 2)

# ms style map -----
# Step 1: Check original data
#print("=== STEP 1: CHECKING ORIGINAL DATA ===")
#print(head(plot_data))
#print(summary(plot_data[c("x", "y")]))
#print(unique(plot_data$category))

# Step 2: Transform to sf object
congruence_sf <- st_as_sf(
  plot_data, 
  coords = c("x", "y"), 
  crs = mcbryde_thomas_2  # Assuming WGS84
)

# Step 3: Project to McBryde-Thomas 2
congruence_sf <- st_transform(congruence_sf, crs = mcbryde_thomas_2)

# Step 4: Check world projection
world_projected_congruence <- st_transform(world, crs = mcbryde_thomas_2)

# Step 5: Create globe border
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))
globe_border_congruence <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Step 6: Check intersection
points_in_globe <- st_intersects(congruence_sf, globe_border_congruence, sparse = FALSE)

# Step 7: Extract coordinates for verification
coords <- st_coordinates(congruence_sf)

# Step 8: Create theme
my_theme_congruence <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank(),
    axis.text = element_blank(),      # Add this line
    axis.ticks = element_blank()      # Add this line
  )

# Step 9: Create plot with debugging
congruence_map <- ggplot() +
  geom_sf(data = congruence_sf, aes(color = category), size = 1.2, alpha = 0.85) +
  geom_sf(data = world_projected_congruence, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_congruence, fill = NA, color = "lightgrey", size = 0.5) +
  scale_color_manual(
  values = c("All three indices" = "#8B0000",
            "IUCN + FUSE" = "#FFD700",           # Changed from "FUSE + IUCN"
            "IUCN + EDGE2" = "#DC143C",          # Changed from "EDGE2 + IUCN"
            "FUSE + EDGE2" = "#FF8C00",          # Changed from "EDGE2 + FUSE"
            "IUCN only" = "#9370DB",
            "FUSE only" = "#32CD32",
            "EDGE2 only" = "#4169E1"),
  breaks = c("All three indices",
             "IUCN + FUSE", 
             "IUCN + EDGE2", 
             "FUSE + EDGE2",
             "IUCN only", 
             "FUSE only", 
             "EDGE2 only"),              # Finally: FUSE + EDGE2
    name = "High priority congruence",
    guide = guide_legend(
      override.aes = list(size = 4, alpha = 1),
      title.position = "top", 
      title.hjust = 0.5,
      ncol = 4
    )) +
  labs(
    title = "Global conservation priority congruence analysis",
    subtitle = "High priority areas (>0.9) showing agreement patterns across IUCN, FUSE and EDGE2 prioritisations",
    x = NULL, 
    y = NULL
  ) +
  my_theme_congruence

ggsave(here::here("outputs", "congruent_high_prioritiy_areas_09_all_indices.png"), 
       congruence_map, 
       width = 10, height = 8, dpi = 300, bg = "white")


# Display both plots
print(congruence_map)

Congruence tests

# Statistical congruence analysis following the methods section
# First, extract priority values for the comparison
priority_threshold <- 0.9

# Create binary matrices for high-priority areas (>0.9)
fuse_high <- CAPTAIN2_FUSE_data_nonzero$Priority > priority_threshold
edge2_high <- CAPTAIN2_EDGE2_data_nonzero$Priority > priority_threshold
iucn_high <- CAPTAIN2_IUCN_data_nonzero$Priority > priority_threshold

# Calculate observed overlaps
observed_fuse_edge2 <- sum(fuse_high & edge2_high)
observed_fuse_iucn <- sum(fuse_high & iucn_high)
observed_edge2_iucn <- sum(edge2_high & iucn_high)
observed_all_three <- sum(fuse_high & edge2_high & iucn_high)

# Calculate minimum high-priority cells for each comparison
min_fuse_edge2 <- min(sum(fuse_high), sum(edge2_high))
min_fuse_iucn <- min(sum(fuse_high), sum(iucn_high))
min_edge2_iucn <- min(sum(edge2_high), sum(iucn_high))

# Calculate percentage overlaps (conservative estimate)
overlap_fuse_edge2 <- (observed_fuse_edge2 / min_fuse_edge2) * 100
overlap_fuse_iucn <- (observed_fuse_iucn / min_fuse_iucn) * 100
overlap_edge2_iucn <- (observed_edge2_iucn / min_edge2_iucn) * 100

# Randomization test function
randomization_test <- function(index1_high, index2_high, n_permutations = 999) {
  observed_overlap <- sum(index1_high & index2_high)
  min_cells <- min(sum(index1_high), sum(index2_high))
  observed_percentage <- (observed_overlap / min_cells) * 100
  
  # Generate random overlaps
  random_overlaps <- replicate(n_permutations, {
    # Randomly shuffle the spatial distribution of index1
    shuffled_index1 <- sample(index1_high)
    random_overlap <- sum(shuffled_index1 & index2_high)
    (random_overlap / min_cells) * 100
  })
  
  # Calculate p-value
  if (observed_percentage > mean(random_overlaps)) {
    p_value <- sum(random_overlaps >= observed_percentage) / n_permutations
  } else {
    p_value <- sum(random_overlaps <= observed_percentage) / n_permutations
  }
  
  return(list(
    observed_percentage = observed_percentage,
    expected_percentage = mean(random_overlaps),
    p_value = p_value,
    random_overlaps = random_overlaps
  ))
}

# Perform randomization tests
print("=== STATISTICAL CONGRUENCE ANALYSIS ===")

[1] "=== STATISTICAL CONGRUENCE ANALYSIS ==="

cat("Threshold for high-priority areas:", priority_threshold, "\n\n")

Threshold for high-priority areas: 0.9 

# FUSE vs EDGE2
test_fuse_edge2 <- randomization_test(fuse_high, edge2_high)

# FUSE vs IUCN
test_fuse_iucn <- randomization_test(fuse_high, iucn_high)

# EDGE2 vs IUCN
test_edge2_iucn <- randomization_test(edge2_high, iucn_high)

# Create a summary table
congruence_results <- data.frame(
  Comparison = c("FUSE vs EDGE2", "FUSE vs IUCN", "EDGE2 vs IUCN"),
  Observed_Overlap_Percent = c(test_fuse_edge2$observed_percentage,
                               test_fuse_iucn$observed_percentage,
                               test_edge2_iucn$observed_percentage),
  Expected_Overlap_Percent = c(test_fuse_edge2$expected_percentage,
                               test_fuse_iucn$expected_percentage,
                               test_edge2_iucn$expected_percentage),
  P_Value = c(test_fuse_edge2$p_value,
              test_fuse_iucn$p_value,
              test_edge2_iucn$p_value),
  Significance = c(
    ifelse(test_fuse_edge2$p_value < 0.001, "***",
           ifelse(test_fuse_edge2$p_value < 0.01, "**",
                  ifelse(test_fuse_edge2$p_value < 0.05, "*", "ns"))),
    ifelse(test_fuse_iucn$p_value < 0.001, "***",
           ifelse(test_fuse_iucn$p_value < 0.01, "**",
                  ifelse(test_fuse_iucn$p_value < 0.05, "*", "ns"))),
    ifelse(test_edge2_iucn$p_value < 0.001, "***",
           ifelse(test_edge2_iucn$p_value < 0.01, "**",
                  ifelse(test_edge2_iucn$p_value < 0.05, "*", "ns")))
  )
)

print(congruence_results)

     Comparison Observed_Overlap_Percent Expected_Overlap_Percent    P_Value
1 FUSE vs EDGE2                 94.56324                 94.89692 0.05405405
2  FUSE vs IUCN                 74.34031                 75.60502 0.00000000
3 EDGE2 vs IUCN                 95.18378                 94.83099 0.01101101
  Significance
1           ns
2          ***
3            *

Relationship with fishing effort : bivariate maps

# Load your CAPTAIN data
IUCN_data <- CAPTAIN2_IUCN_data_nonzero  # Assuming this exists from your previous code
FUSE_data <- CAPTAIN2_FUSE_data_nonzero  # Assuming this exists from your previous code
EDGE2_data <- CAPTAIN2_EDGE2_data_nonzero # Assuming this exists from your previous code

# Process and prepare the data function
process_bivariate_data <- function(priority_data, fishing_data, min_priority = 0) {
  # Filter for cells with priority > min_priority
  priority_data_filtered <- priority_data %>%
    filter(Priority > min_priority)
  
  # Prepare fishing data
  fishing_data <- fishing_data %>%
    rename(Longitude = lon_05deg, Latitude = lat_05deg, FishingHours = predicted_fishing_hours)
  
  # Join datasets
  combined_data <- priority_data_filtered %>%
    left_join(fishing_data, by = c("Longitude", "Latitude"))
  
  # Handle NAs in fishing hours (replace with 0)
  combined_data$FishingHours[is.na(combined_data$FishingHours)] <- 0
  
  # Normalize priorities to 0-1 range and log transform fishing hours
  combined_data <- combined_data %>%
    mutate(
      Priority_Norm = Priority, # Assuming already in 0-1 range
      # Log transform fishing hours to better handle skewed distribution
      FishingHours_Log = log1p(FishingHours),
      # Normalize log fishing hours
      FishingHours_Norm = (FishingHours_Log - min(FishingHours_Log, na.rm = TRUE)) /
                         (max(FishingHours_Log, na.rm = TRUE) - min(FishingHours_Log, na.rm = TRUE))
    )
  
  # Set projection
  mcbryde_thomas_2 <- "+proj=mbt_s"
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process the three datasets with fishing data (only cells with Priority > 0)
iucn_fishing_sf <- process_bivariate_data(IUCN_data, aggregated_data, min_priority = 0)
fuse_fishing_sf <- process_bivariate_data(FUSE_data, aggregated_data, min_priority = 0)
edge2_fishing_sf <- process_bivariate_data(EDGE2_data, aggregated_data, min_priority = 0)

# Create color palette - using the PurpleOr scheme
create_bivariate_palette <- function() {
  map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
  map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
  map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
  map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
  map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
  map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
  map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
  map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
  map_pal_mtx[1, 1] <- '#ffffee'
  map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))
  return(map_pal)
}

map_pal <- create_bivariate_palette()

# Color mapping function
get_color <- function(priority, fishing) {
  priority_class <- cut(priority, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fishing_class <- cut(fishing, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fishing_class)-1)*4 + as.numeric(priority_class)])
}

# Apply colors to the datasets
iucn_fishing_sf$bivariate_color <- mapply(get_color,
                                          iucn_fishing_sf$Priority_Norm,
                                          iucn_fishing_sf$FishingHours_Norm)

fuse_fishing_sf$bivariate_color <- mapply(get_color,
                                          fuse_fishing_sf$Priority_Norm,
                                          fuse_fishing_sf$FishingHours_Norm)

edge2_fishing_sf$bivariate_color <- mapply(get_color,
                                          edge2_fishing_sf$Priority_Norm,
                                          edge2_fishing_sf$FishingHours_Norm)

# Get world map and project
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create globe border
lon_points <- rep(c(-180, 180), each = 100)
lat_points <- c(seq(-90, 90, length.out = 100), seq(90, -90, length.out = 100))
globe_outline <- data.frame(lon = lon_points, lat = lat_points) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_combine() %>%
  st_cast("POLYGON") %>%
  st_transform(crs = mcbryde_thomas_2)

# Function to create individual bivariate maps with labels and titles
create_bivariate_map <- function(data_sf, label, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = bivariate_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_outline, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL, title = title) +  # Added title
    theme(panel.grid = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(size = 14, face = "bold", hjust = 0.5),  # Title styling
          plot.margin = margin(5, 5, 5, 5, "mm")) +
    # Add smaller label in top-left corner
    annotate("text", x = -Inf, y = Inf, label = label,
             hjust = -0.2, vjust = 1.2, size = 4, fontface = "bold")  # Reduced from 8 to 4
}

# Create the three maps with labels and titles
iucn_plot <- create_bivariate_map(iucn_fishing_sf, "(A)", "IUCN")
fuse_plot <- create_bivariate_map(fuse_fishing_sf, "(B)", "FUSE")  
edge2_plot <- create_bivariate_map(edge2_fishing_sf, "(C)", "EDGE2")

# Create a larger legend for better visibility
larger_legend <- bi_legend(pal = map_pal, dim = 4,
                          xlab = 'Conservation\npriority',
                          ylab = 'Fishing\neffort',
                          size = 4) +  # Increased from 3 to 4 for bigger legend
  theme(
    axis.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    axis.text = element_blank(),
    legend.text = element_text(size = 8),  # Reduced from 10 to 8
    legend.title = element_text(size = 10, face = "bold"),  # Reduced from 12 to 10
    plot.margin = margin(10, 10, 10, 10, "mm")
  )

# Create layout matrix: 3 rows, 2 columns
# Column 1: maps (wider)
# Column 2: empty, legend (middle row), empty
layout_matrix <- rbind(
  c(1, NA),      # Row 1: Map A, empty
  c(2, 4),       # Row 2: Map B, legend
  c(3, NA)       # Row 3: Map C, empty
)

# Arrange all plots together (no auto-display)
combined_plot <- gridExtra::arrangeGrob(
  iucn_plot,
  fuse_plot,
  edge2_plot,
  larger_legend,
  layout_matrix = layout_matrix,
  heights = c(1, 1, 1),      # Equal height for all three map rows
  widths = c(3, 1)           # Maps take 3/4 width, legend takes 1/4
)

# Save with dimensions better suited for A4 portrait
ggsave(here::here("Outputs", "Fig.4_All_indices_fishing_bivariate_maps_vertical.png"),
       combined_plot, width = 8.3, height = 11, dpi = 300, bg = "white") 

# Display the saved image in the Quarto document
knitr::include_graphics(here::here("Outputs", "Fig.4_All_indices_fishing_bivariate_maps_vertical.png"))